Bose-Einstein condensation in trapped bosons: A Variational Monte Carlo analysis 
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Several properties of trapped hard sphere bosons are evaluated using variational Monte Carlo 
techniques. A trial wave function composed of a renormalized single particle Gaussian and a hard 
sphere Jastrow function for pair correlations is used to study the sensitivity of condensate and non- 
condensate properties to the hard sphere radius and the number of particles. Special attention is 
given to diagonalizing the one body density matrix and obtaining the corresponding single particle 
natural orbitals and their occupation numbers for the system. The condensate wave function and 
condensate fraction are then obtained from the single particle orbital with highest occupation. 
The effect of interaction on other quantities such as the ground state energy, the mean radial 
displacement, and the momentum distribution are calculated as well. Results are compared with 
Mean Field theory in the dilute limit. 



I. INTRODUCTION 

The spectacular demonstration of Bose-Einstein con- 
densation (BEC) in gases of alkali atoms '^^Na, '^Li 
confined in magnetic traps ^, ^ has led to an explosion 
of interest in confined Bose systems. Of interest is the 
fraction of condensed atoms, the nature of the conden- 
sate, the excitations above the condensate, the atomic 
density in the trap as a function of Temperature and the 
critical temperature of BEC, Tc- The extensive progress 
made up to early 1999 is reviewed by Dalfovo et al.p. 

A key feature of the trapped alkali and atomic hydro- 
gen systems is that they are dilute. The characteristic di- 
mensions of a typical trap for ^'^Rb is Uho — {fi/muj^)^ — 
1 - 2 X 10^ A (Ref. 1). The interaction between Rb 
atoms can be well represented by its s-wave scattering 
length, a Jib- This scattering length lies in the range 
85 < am < 140ao where ao = 0.5292 A is the Bohr 
radius Q. The definite value a^b = lOOao is usually se- 
lected and for calculations the definite ratio of atom size 
to trap size aRb/am ~ 4.33 x 10~^ is usually chosen Gj. 
A typical ^"^ Rb atom density in the trap is n ~ 10^^ — kF"* 
atoms/cm'^ giving an inter-atom spacing £ ~ 10^ A. Thus 
the effective atom size is small compared to both the trap 
size and the inter-atom spacing, the condition for dilute- 
ness (i.e., na|jj ~ 10~^ where n = N/V is the number 
density). In this limit, although the interaction is impor- 
tant, dilute gas approximations such as the Bogoliubov 
theory Q, valid for small na^ and large condensate frac- 
tion no = Nq/N, describe the system well. Also, since 
most of the atoms are in the condensate (except near Tc), 
the Cross-Pitaevskii equation[^ |^ for the condensate de- 
scribes the whole gas well. Effects of atoms excited above 
the condensate have been incorporated within the Popov 
approximation^. One of the chief purposes of this pa- 
per is to go beyond the dilute limit, to test the limits 
of these approximations and to explore the properties of 
the trapped Bose gas as na"^ increases between the di- 
lute limit and the dense limit. We use Variational Monte 
Carlo (VMC) methods. We increase the density by in- 
creasing both N and the s- wave scattering length up to 
the value na^ ~ 0.21 which describes liquid ^He at SVP 



when the ^He atoms are represented by hard spheres of 
diameter a = 2.203 A||l|l. 

In addition to the mean-field theories noted above, the 
trapped Bose gas at finite temperatures has been inves- 
tigated using Path Integral Monte Carlo (PIMC) meth- 
ods. KrauthQ simulated 10,000 atoms in a spherical trap 
with ratio of scattering length a to trap length given by 
a/o/io = 4.3 X 10"'^ noted above. He showed that the crit- 
ical temperature Tc is lowered compared to the ideal Bose 
gas as a result of interaction. Tc is lower because the re- 
pulsion between the atoms spreads the atoms in the trap 
and lowers the density compared to the non- interacting 
case. The same result has been obtained in mean-field 
approximations Krauth also showed that, in the di- 
lute limit, the condensed atoms are highly concentrated 
at the center of the trap while the uncondensed or ther- 
mal atoms are spread out over a wide range, are dilute 
and well approximated by a classical, ideal gas. There 
is little interaction between the condensed and uncon- 
densed components (see also 0). 

Criiter et al.[^ evaluated Tc for a uniform, bulk, hard 
sphere Bose gas over a wide density range, from dilute 
to hquid ^He densities. They find that Tc is increased 
above the ideal Bose gas value by interaction in the dilute 
range. In the uniform gas case the density is not changed 
by interactions. At liquid ^He densities, Tc is decreased 



by interaction [p, 7 
Holzmann et a 



. ||ll]| made a direct comparison between 
Hartree-Fock (HP) and PIMC calculations of the number 
density N{r) of atoms in a trap. The atoms were again 
represented by hard spheres with a/a^o — 0.0043. From 
N(r), they find for temperatures near Tc that the conden- 
sate fraction Nq is larger in the exact PIMC evaluation 
than in the HP approximation. The energy beyond the 
Hartree-Fock approximation is often denoted the "corre- 
lation" energy. This correlation apparently allows iVo to 
increase at a given T and allows condensation to begin at 
a higher temperature. At lower temperature T < 0.75Tc, 
there is excellent agreement between the PIMC and HP 
N{r). The increase in Nq with exact representation of 
the interaction effects is consistent with the correspond- 
ing increase in Tc with interaction in the uniform Bose 
gas. 
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Giorgini et al. [|2j have evaluated the ground state en- 
ergy E/N and the condensate fraction Nq/N at T = K 
of a uniform Bose gas over a wide density range (10^^ < 
na^ < 10-1) using Green Fimction Monte Carlo (GFMC) 
methods. They find that the mean field results of E/N 
and the Bogoliubov result for Nq/N agree well with the 
GFMC values in the density range 10"^ < na^ < 10~^. 
However, there are clear differences at higher densities 
na^ > 10~^ (helium density is na^ ~ 0.21). The re- 
sults are not sensitive to reasonable variation of the inter- 
boson potential p^. 

In this context, we have evaluated the ground state 
properties of a trapped, hard sphere Bose gas over a 
wide range of densities using Variational Monte Carlo 
(VMC) methods. We begin in the dilute limit, small N 
and a/atio = 0.0043 corresponding to ^"^ Rh in a trap, 
and increase both N and a separately to increase the 
density up to liquid ^He densities. At the lower densi- 
ties we compare the energy E/N and root mean square 



amplitudes < x 



V 



^ > and < 2;^ > of atoms in an 



anisotropic trap with Gross-Pitaevskii (GP) results [|18[. 
The two methods agree well at low densities but even 
at densities na^ « lO^'^ small differences in E/N are 
readily apparent. Also, at higher densities we find the 
effects of interaction depend separately on N and a/uho, 
not simply in the product Na/uho as appears in the GP 
theory. As density is increased still further, we find that 
the condensate is no longer concentrated at the center of 
the trap. Rather increased interaction (increased a/aho) 
depletes the condensate at the center and the conden- 
sate appears at the edges of the trap - as found in liq- 
uid ''He droplets (Lewart et al.[|l3j). Also, as density in- 
creases, correlations in the single particle density appear, 
reflecting the interaction in a confined space, in the same 
way that interaction at higher density introduces correla- 
tion in the pair correlation function in the uniform case. 
We evaluate the momentum distribution and condensate 
fraction over the dilute to dense range and compare with 
mean- field results p9[|. 

In section 2 we introduce the Hamiltonian, the wave 
function and the MC method and the definition of the 
natural orbitals. The results are presented in section 3 
and discussed in section 4. 



II. BOSONS IN A HARMONIC TRAP 



We consider both a spherically symmetric (S) harmonic 
trap and an elliptical (E) harmonic trap. 
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K.t(r) = 



iS) 
(E) 



(2) 



Here usf^^ defines the trap potential strength. In the case 
of the elliptical trap, Vexti^^y, z), uj^o = is the trap 
frequency in the perpendicular or xy plane and lOz the 
frequency in the z direction. The mean square vibra- 
tional amplitude of a single boson at T = QK in the trap 
(^) is < x^ >= (h/2mu!ho) so that aho = {'h/mLOho)'^ 
defines the characteristic length of the trap. The ratio of 
the frequencies is denoted A = w^/w^ leading to a ratio 
of the trap lengths (aj_/az) = {i^z/'^i.)^ = \/A. 

We represent the inter boson interaction by a pairwise, 
hard core potential 



oo r < a 
r > a 



(3) 



where a is the hard core diameter of the bosons. Clearly, 
Vint (f) is zero if the bosons are separated by a distance r 
greater than a but infinite if they attempt to come within 
a distance r < a. 

The weak interaction limit is a << a^o and a << n~3 
(where n = N/V is the local number density), a hard core 
diameter small compared to the dimensions of the trap 
and compared to the inter-particle spacing I = {V/N)^ . 
For trapped alkali atoms we have typically na?< 10"^. 
Introducing lengths in units of a/jo, r — > r / a^o, and fiLOho 
as units of energy as in Q, the Hamiltonian is: 

^ 1 

^ = E 2(-v'+^'+2/'+A'^,')+E ^»*(k.-r,l). (4) 

i i<j 

Since there is Bose condensation, we have nXj^> 2.616 
where At is the atomic thermal wavelength. Thus we are 
in the regime where the atomic wavelength is long com- 
pared to the hard core diameter, Xt » a or ka << 1 
where k = I-k/Xt- The scattering of two particles in- 
teracting via a hard core potential in the limit fca << 1 
is purely s wave with scattering length a. If we approx- 
imate the full potential between the two particles by a 
contact potential. 



A. The System 



v(r) = g5(y) 



Sir) 



(5) 



We consider N bosons of mass m confined in an ex- 
ternal trapping potential, Vext{i^), and interacting via a 
two-body potential Vint (ri ,^2)- The Hamiltonian for this 
system is: 



N 



^ = E ( + Vextir^) ) + J2 r,). (1) 



N 



the scattering length in this limit between the two is 
again purely s wave with scattering length a. Thus we 
may compare directly results calculated using a hard core 
potential (3) and with a contact potential approximation 
(5) in the regime a « Xt- Specifically, we may com- 
pare the present MC results with results calculated using 
(5) and the mean field. Gross- Pitaevskii (GP) equation. 
This comparison is especially interesting in the dilute 
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limit na^ << 1. At high densities {na^> 0.1), we ex- 
pect short range, pair correlations induced by the hard 
core to be important and short range correlations are not 
well described by a mean field theory. 



B. Wave Function 

To describe the ground state of the N bosons, we intro- 
duce a variational trial wave function which is a product 
of a single particle function g{r) and a pair Jastrow func- 
tion II /(In -ral), 

^'^(r,..r„,a,/3) ^ Y[g{a, /3,r,)'[[f{a,\r, - r^l) (6) 

i i<j 

where a and (3 are the variational parameters. We select 
a single particle function, 

g{a, P, r,) = exp [~a{x'^ + yf + (3x1)] (7) 

which is a HO ground state function having two varia- 
tional parameters, a and /?. For spherical traps, [3—1, 
and for noninteracting bosons (a = 0), a = \/2a\^. 
For the pair function we select the exact solution of the 
Schroedinger equation for two particles interacting via 
the hard core potential (3) in the limit fc ^ 0, i.e. 

(see for example Huang (2|]). The 5',y(rj^ ..r„) therefore 
has the correct form for small jr^ — Yj\ and has two vari- 
ational parameters a and [3 that describe the spread of 
the bosons in the trap as the hard core diameter is in- 
creased. By constructing the wave function in this way, 
we limit the number of variational parameters while pre- 
serving the correct functional form in the a — > limit. 
However, the lack of any variational parameters in the 
Jastrow term is a potential source of inaccuracy. 

We then minimize the expectation value of the Hamil- 
tonian as obtained from 

with respect to a and /? using the Metropolis Monte Carlo 
method of integration. This is accomplished by using a 
Metropolis random walk to generate a set of M par- 
ticle configurations, fli..ftM which conform to the prob- 
ability distribution l^'i/p. We then approximate < H > 
by summing over the 'local energy' as follows 

n—l ^ \ J / 

For a review of the Variational Monte Carlo method, see 




C. Condensate and Natural Orbitals 

A goal here is to calculate the condensate fraction and 
condensate density in the ground state. To do this we 
require a definition of the condensate single particle state. 
Following Penrose and Onsager, Lowdin and others |2^, 
wc take the one-body density matrix (OBDM) as 
the fundamental quantity for an interacting system and 
define the natural single particle orbitals (NO) in terms 
of the OBDM. 

The OBDM is ^ 

p{y',y)^<¥{y'),^{y)>, (9) 

where 4'(r) is the field operator that annihilates a single 
particle at the point r in the system. At T = OK, the 
expectation value is evaluated using the wave function 
in (6). To define the NO, we introduce a set of 
single particle states having wave functions (/>i(r) and ex- 
pand ^(r) in terms of the operators di which annihilate 
a particle from state |i >, 

*o = ^'^»(r)a». (10) 

i 

Requiring that the satisfy the usual commutation 
— and number relations (< a\aj >— NiSij), 

we have 

p{r,Y')^J:^,r,{r')Mr)NAJ 

This may be taken as the defining relation of the NO, 
(j)i{Y). Specifically, we have from (11), 

J dY..dY'4>*{Y)p{Y, y')c^,{y') = NAj, (12) 

SO that the NO may be obtained by diagonalizing the 
OBDM. The eigenvectors are the NO and the eigenvalues 
are the occupation, Ni, of the orbitals. The condensate is 
the orbital having the highest occupation, denoted 0o ('') 
and the condensate fraction is no = Nq/N. 

The relations (11) and (12) involve the vector r and r' 
and cannot be solved directly as matrix equations. To 
obtain matrix equations, we restrict ourselves to spher- 
ical traps and seek equations for the radial component 
of the NO. Assuming the potential seen by a single par- 
ticle, including inter-particle interaction, is spherically 
symmetric, the NO will have the form 

UY)^Rniir)Yimie,(b) (13) 

Where Y/m(f^) are the spherical harmonics and Rni{r) 
is the radial wave function and i = n,l,m are the state 
indices. We expand the OBDM in its angular momentum 
(l) components pi(r,r') as 

p(r,r') = Yl ^^^Pi(^ ■ i')Piiry) (14) 
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where Pi{cos9) are Legendre Polynomials in the angle 
between f and f' and 



(15) 



Substituting (13) and (14) into (11) and using the prop- 
erties of the spherical harmonics, Yim{^)^ we obtain a 
relation equivalent to (11) for each I component 4>ni{r) 
of the NO, 



(16) 



where i = nl. To solve this equation readily as a matrix 
equation, we introduce the radial function 

Uni{r) = r<j)ni{r). (17) 

The Uni{r) is well behaved at r ^ and has dimensions 
like a one-dimensional wave function. In terms of 
the Uni{r), the defining relations are: 



[rp, {r, , r[ )r'] = ^ u^i (r)<, (r') A^„i 



(18) 



and 



dr[rp,{r,,r'y]u„i{r') - Uni{r)Nni (19) 



The (r^, r' )r'] serves as a one-dimensional OBDM 
along the radial coordinate. This ID matrix relation may 
be solved numerically on a grid in r to obtain the Uni (r) 
as eigenvectors and the Nni as eigenvalues. The conden- 
sate orbital is 

Mr) - -^uooir)/r. (20) 



Air 

The momentum distribution may be obtained from 
the OBDM as well given that the orbitals in momentum 
space are 

7 „ . 1 



dre 



and the momentum distribution is 

p(k)=^n,|0,(k)|2. 



(21) 



(22) 



Substituting ( pl| ) into ( |2^ ) obtains an expression for /5(k) 
in terms of the OBDM: 



P(^^ = J^ J drdr'p{r,r')e 



/^p^k•(r-r') 



(23) 



The momentum distribution for a spherically symmetric 
system is therefore 

m = m) = 

i /rfri-rfr. mr^-r^r ^^4) 
X f°dr ^i^(i'i+r, r2-.rAr) 

where the orientation of r may be chosen arbitrarily. 
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FIG. 1: The condensate orbital, (po{r), obtained by numerical 
diagonalization of the OBDM po(r,r') = ^ - riit/i* (r)0i (r') calcu- 
lated by Variational Monte Carlo (VMC) (solid dots) for an ideal 
Bose gas in a harmonic trap. The dashed line is the harmonic os- 
cillator (HO) ground state wave function for the same trap. 4>o{r) 
and r are dimensionloss in units of a/i^. 



III. RESULTS 

In this section we present the results of the present 
Variational Monte Carlo (VMC) calculation of proper- 
ties of bosons at T = K confined in a harmonic trap. 
To begin, Fig. 1 shows the condensate orbital (t>o{r) for 
independent, non- interacting bosons in a spherical har- 
monic trap compared to the harmonic oscillator (HO) 
ground state function for the same trap. The (l)o{r) is 
calculated by evaluating the one body density matrix 
(OBDM) and diagonalizing it numerically to obtain the 
single particle orbitals and their occupation as discussed 
in section 2. For no interaction, only 0o(^) is occupied 
(?io = 1)- The excellent agreement between 0o('') and the 
HO ground state function for all r provides a good check 
of the method for non-interacting bosons. The statistical 
sample is proportional to and the statistics become 
poor at r — !■ 0. 

Fig. 2 shows the energy per particle, E/N, of a gas of 
N weakly interacting hard sphere bosons in an ellipsoidal 
trap. The ratio of the characteristic length of the short 
axis (z axis) to the longer perpendicular [x and y) axis of 
the trap is (az/zj_) = I/VA;, A = VS- The hard sphere 
diameter, a, (scattering length) of the bosons corresponds 
to atoms in an ellipsoidal trap with a/aho = 0.00433 
with a± = Uho- In Fig- 2 the dots are the present 
VMC E/N of the whole Bose gas while the dashed line 
is the E/N oi the condensate calculated by Dalfovo and 
Stringari||l^ using the Gross-Pitaevskii (CP) equation. 
Our purpose is to make comparisons between the present 
VMC calculation and CP equation results across the di- 
lute regime for which CP is expected to be valid. The 
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FIG. 2: Energy per particle, in units of huij^, for hard sphere 
bosons in an anisotropic trap as a function of the number of par- 
ticles N in the trap. Solid circles are the present VMC values for 
hard spheres with diameter corresponding to the scattering length 
of Rb, anij = 4.33 X lO^^a^Q - where a^Q is the trap length in the 
perpendicular direction. The error bars lie within the solid dots. 
The dashed line is the value obtained using the Gross-Pitaevskii 
equation (GP) for the same system [h 8| . 



region 100 < < 20000 corresponds to atom densities 
at the center of the trap of 2 x 10"^ < na^ < 2 x IQ-^. 
In the earhest experiments with Rb there were typ- 
ically N = 10, 000 atoms in the trap. In more recent 
experiments N is larger, TV ~ 10^ — 10^. 

At small N values in the dilute limit, there is excel- 
lent agreement between the GP and VMC energies. In 
this regime, the mean-field, GP equation is expected to 
be accurate. As N increases a clear difference between 
the VMC and GP E/N values emerges. We find that 
for a fixed scattering length, the difference 6{E/N) = 
[Emc ~ Ecp)/N is proportional to A^s and can be well 
represented as 

5{E/N) = [(3.0xlO"^±10"^)x(iV)t-0.052]?kj_L. (25) 

The difference is 1.8% at iV = 10,000 and 2.5% at 
N — 20, 000. The difference, we believe, arises largely 
because there is excitation of bosons above the conden- 
sate. Emc/N includes the excited atoms while Egp/N 
is the energy of the condensate alone. We return to this 
point in the Discussion. 

Fig. 3 compares the root mean square displacement 
of hard sphere bosons from the center of the same 
anisotropic trap discussed in Fig. 2 calculated using VMC 
and the GP equation. The upper (lower) line is the radial 
displacement along the _L (z) direction. The agreement 
between the GP and the VMC displacements is excellent, 
right up to iV = 20, 000 for a = am- 
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N 

FIG. 3: Average axial displacement of bosons from the center of 
an anisotropic trap in the ± direction, <x^+y^>'^ , (top) and in 
the z direction <2;^>2 (bottom) expressed in units of the perpen- 
dicular trap length a/^o as a function of the number of particles, 
TV. Solid dots are from the present VMC results for hard spheres 
with diameter am, = 4.33 X 10~^aho ■ The dashed lines are the 
values obtained from the Gross-Pitaevskii equation (GP) for the 
same a^j, |18|. 



It is interesting that the VMC and GP displacements 
agree well while the energies in Fig. 2 differ for N> lO'*. 
Essentially, the E/N is very sensitive to the few high 
energy bosons at the edges of the trap. That is, the small 
number of atoms having large displacements increase the 
energy significantly but change < > little. An example 
of this effect is found in the Thomas- Fermi approximation 
where the density is cut off to zero at a specific radius [|| . 
In the TF model there is no tail in the density reaching 
up to large r values. The cut off changes < > little 
but E/N is significantly affected Q. 

Fig. 4 again shows E/N for hard sphere bosons in an 
anisotropic trap as a function of Na/aho- However, in 
this case the product Na/auo is adjusted by varying both 
N and a/auo- The star symbols show E/N for a = aju, 
and 2000 < < 20,000, the crosses for a = lOa^b 
and 200 < < 2, 000 and the square a = 20am, and 
100 < A^ < 1, 000. If the impact of interaction and E/N 
depended solely on the product Na/aho, as is the case 
in the mean-field, GP equation, all three lines in Fig. 4 
would coincide. E/N clearly depends separately on A^ 
and a/ttfio, even in the region of A^ = 10,000 — 20,000. 
The separate dependence is not large at Na/aho ~ 20 but 
becomes increasingly large as Na/aho increases. Also, at 
these and larger densities, the parameter which deter- 
mines the magnitude of corrections arising from interac- 
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FIG. 4: Energy per particle, in units of hui^, for interacting bosons 
in an anisotropic trap as a function of the product Na/a/ic, (number 
of particles N and scattering length a) as obtained by the present 
VMC calculation for hard spheres with core diameter a = 1, 10, 20 
times the scattering length of *^Rb, aju, = 4.33 X 10~^af^o- Open 
squares are Gross-Pitaevskii equation results for the same Na. Es- 
timated errors lie within symbol size. 
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FIG. 5: Condensate natural orbital <^o('') for 128 hard spheres 
with core diameter a = 0,1, ...,64 times the scattering length of 
^'^Rb, apii, = 4.33 X 10~'^ahQ. in a spherically symmetric harmonic 
trap All lengths are in terms of the trap length a^Q = h/muif^o- 



tions is N 5 {a / ttho)'^ ■ Apparently the interaction effects 
depending on Na/aho are valid only in the limit of small 
a {a/afio << 1). We examine the functional form of the 
separate dependence of E/N on A'^ and a/a^o in the dis- 
cussion section. 

Having investigated lower densities and made com- 
parisons with results obtained using the GP equation, 
we now turn to higher densities and bosons represented 
by hard spheres having larger hard core diameters. We 
evaluate the OBDM and the density for these cases go- 
ing up to densities comparable to liquid *He droplets. 
Fig. 5 shows the condensate orbital (wave function) 
for 128 bosons in a spherical harmonic trap as the 
hard core radius, a, is increased from zero up to a = 
64ai^b {aRb/cLho = 0.00433). The case a/ ant = 1 corre- 
sponds to A'' = 128 Rb atoms in a spherical trap. Clearly 
the condensate orbital spreads out in the trap as a in- 
creases. At a/afO) = 64, the condensate density is effec- 
tively constant in the trap out to nearly three times the 
trap length parameter aho- For these larger core radii, 
the appropriate measure of the interaction is na^ where 
n = N/V. For a — 64afli,, na^ ~ 2 x 10^^ at the center 
of the trap. 

In Fig. 6 we show the total density p{r) and the 
density of atoms in the condensate orbital as A^|(/)o(r)p 
for 64 bosons in a spherical trap with a increased to 
128afli,, 256afli, and 512afli,. Since we plot A'|(/)o(r)p 
rather than Ao|(/)o(^)Pj "condensate density" can exceed 
the total density. At a = 256afli, the condensate is 
'^0 = Nq/N = 20%. We note that as a increases, the 
condensate moves away from the center of the trap. At 
a = 512ai{h, the condensate is at the edges of the trap. 
When na^ is large at the center of the trap, the maxi- 



mum condensate density is in the region of lower parti- 
cle density found at the edge of the trap, as calculated 
for liquid "'He droplets Thus the location of the 

condensate is entirely different at small and large scat- 
tering length. In a slave boson approach, depletion of 
the condensate at the center of the trap has also been 
demonstrated j2^. Also, at large a /aho the total den- 
sity develops correlations. These correlations reflect the 
inter-boson correlations induced by the hard core inter- 
action. For a = 512a/fft we have a/aho — 2.2. With a 
peak in the density at r = 0, we expect the first mini- 
mum in the density at r ~ a/aho — 2.2 as in the upper 
frame of Fig. 6. 

Fig. 7 gives the fraction of bosons in the conden- 
sate orbital, no, calculated by VMC and diagonaliza- 
tion of the OBDM corresponding to the condensate or- 
bitals shown in Fig. 5. The no values, as in Fig. 5, 
are for 128 hard sphere bosons in a spherical trap with 
hard sphere radius a/a^h — 1,2,4,16,32 and 64 e.g. 
a = 64 X Oflb = 64 X 4.33 x 10"^ = 0.277a/jo. The Bo- 
goliubov Q result for no adapted to a spherical trap|T9| 
is shown as a dashed line. Visible departures of the VMC 
no from the Bogoliubov result begin at no ~ 0.96 (a de- 
pletion of 4%) corresponding to a density at the center 
of the trap na^ ~ 1 x 10"'^. In the uniform gas case, the 
Bogoliubov result remains accurate up to a condensate 
fraction no ~ 0.89 (11% depletion) which occurs at a 
density na^ ~ 5 x 10^"^ |l^. The Bogoliubov approxima- 
tion has a more limited range of application for bosons in 
a trap because the interaction changes the density pro- 
file (and the shape of the condensate wave function) as 
well as simply depleting the condensate and no depends 
on the density and shape of the condensate wave func- 
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FIG. 6: Condensate density, |0o('')Pi scaled by the number of 
particles A^, vs. total particle density, N{r), with scattering length 
a = 128, 256, 512 X aju,, from bottom to top — where N{r) is nor- 
malized to N, |i?io('')P is normalized to 1, and all lengths are in 
units of the characteristic trap length ajj^. 



tion. In the uniform case, the density cannot change. In 
Fig 8., we again show the condensate fraction, no, as a 
function of the ratio of the scattering length, a, to the 
characteristic trap length, aho- Here, no is given for three 
different numbers of particles: N = 64, 128 and 256. The 
corresponding value of the particle density, na^, for the 
64 particle case is shown on the top axis for reference. At 
hquid helium densities, the condensate fraction is roughly 
twice that of bulk liquid '^He. This difference can be un- 
derstood by noting the shape of the radial condensate 
density shown for the a = 256aRb case in Fig. 6. Here, 
we see that while the maximum particle density occurs 
at the center the trap, the condensate density is peaked 
in the low density region at the edge of the cloud. This 
dilute region allows for a larger fraction of particles to oc- 
cupy the condensate orbital than in an uniform system 
at "^He densities. 
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FIG. 7: The condensate fraction, no, for 128 particles in a har- 
monic trap as a function of the ratio of scattering length to trap 
length a/a^o using three methods. Occupation numbers for the 
ground state orbital </>o found using Variational Monte-Carlo (MC) 
methods (solid dots) agree well with values of no obtained from the 
Bogoliubov equations for a uniform gas (open ciicles) and a Mean 
Field Bogoliubov approximation (dashed line) [hg| for no > .9. 
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IV. DISCUSSION 



Fig. 9 shows the effect of increased scattering length, a, 
on the momentum distribution of particles in a isotropic 
harmonic trap at T = 0. The values for the scatter- 
ing length and the trap configuration under consideration 
correspond to those shown for the spatial distribution in 
Fig. 5. 



In this section we compare the present MC values for 
the density, condensate fraction and energy of bosons in 
a trap with mean field (MF) and Thomas-Fermi (TF) 
approximation results. The aim is to assess the limits 
of applicability of the MF and TF expressions and to 
investigate the origin of any differences between MF and 
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FIG. 9: Momentum distribution, n(fc), for 128 hard spheres with 
diameter a = 0, 1, ..32 times the scattering length of Rb, a^jj, = 
4.33 X 10~^, in a harmonic trapping potential. Momentum is in 
units of the inverse of the characteristic length of the trap a^J = 



MC values. We start with the density at the center of the 
trap, n(0), since this density is needed in MF expressions 
for the depletion of the condensate and for the energy. 
Also, comparisons [4] of the density calculated in the TF 
approximation and in mean field, Gross-Pitaevskii (GP), 
approximations show that 717^^(0) is accurate for large 
Na/aho {Na/auo ~ 100 ). 

The density of N independent bosons in an asymmetric 
trap is nfto(r) = (^^^(r) = 7V/(7ria:raya^) exp[-(a;2/a^ + 
y^/a^ + z^/a^)]. For the elliptical trap considered above, 

cix — CLy = aho and — aho/V^, the density at the 
center of the trap is [Q 



(26) 



The (j)o{r) for a spherically symmetric trap (A = 1) is 
shown in Fig. 1. As interaction is increased (e.g. a /aho 
is increased) the 4>o{r) spreads out and the density n(0) 
at the center of the trap decreases, as depicted in Fig. 5. 
For A = 1, the ratio of the density at the center with 
interaction to the rihoiQ) for no interaction in the TF 
approximation is 



TTvT = 5 ){Na aho) 

«/io(0) 8 



(27) 



This ratio calculated using MC, nMc(0)/n/io(0), for in- 
creasing scattering length, a/aho = Fa^h = F x 4.33 x 
10~^ with F = 1, 4, 8, 16, 32 and 64 can be obtained from 
Fig. 5 and is listed in Table 1. The TF ratio agrees 
well with the MC ratio in the range 4 < < 16 or 



a/aRb Na/aho 


nMc{0)a^ 


nTF(0)/nfco(0) 


nMc{0)/nho(0) 


1 0.55 


1.3 X 10"'' 


0.93 


0.71 


4 2.22 


4.3 X 10"^ 


0.41 


0.38 


8 4.43 


2.4 X 10"'' 


0.27 


0.25 


16 8.87 


1.0 X 10"^ 


0.17 


0.14 


32 17.7 


4.3 X 10"^ 


0.12 


0.075 


64 35.5 


1.9 X 10"^ 


0.077 


0.041 



TABLE I: The ratio of the density at the center of the trap 
calculated using the Thomas- Fermi expression, nTF(O), and 
in the present VMC evaluation, nMc{0), to the density for 
no interaction, nho(O). The ratio decreases as the scattering 
length / hard core diameter a/ am, is increased as shown in 
Fig. 5. 



2 < Na/aho < 10. We expect the TF result to be inac- 
curate at small Na/aho because the TF limit is a large 
Na/aho approximation. 

In addition, at large Na/aho, the density exceeds the 
dilute limit so that all mean field theories become in- 
accurate. The dilute limit is exceeded at n{0)a^> 10~^ 
which corresponds to > 16 for the gas considered in 
Table 1. Between these limits however, for Na/aho^ 5 
and n{0)a^< 10"'^, the rtTF(O) gives a good estimate of 
the boson number density at the center of the trap. 

For a uniform Bose gas, the MF, Bogoliubov the- 
ory predicts a condensate fraction no = Nq/N = 1 — 
(8/3)(na'^/7r) 2 where n = N/V is the uniform density 
[g. That is, the fraction of bosons excited out of the 
condensate is 



5N _ N -No 



O TT 



(28) 



The corresponding result for bosons in a spherical har- 
monic trap is W 



SN 



57r ,nTF(0)a^ , i 



(29) 



Fig. 7 shows that the mean field expressions (g8|) and (|2£ 
agree well with the present MC values of no for iV = 128 
bosons up to a/aho ~ ^GaRb or up to n{0)a^ « 10""^. 
Specifically, for F = 8 at nMc{0)a^ = 2.4 x lO"-*, the MC 
value for the depletion in Fig. 7 is SN/N — 1.9% while 
the Bogohubov result from (|^) is SN/N = 1.7%. For 
larger values of n{0)a^ the depletion is significantly un- 
derestimated by the Bogoliubov expressions. This agrees 
with the findings of MC determinations for no in a uni- 
form Bose gas||l^. In the uniform case, ( ^8|) agrees with 
MC values up to na'^ sa 5 x 10"'^ and underestimates de- 
pletion at larger na^ . The MF results are accurate up to 
somewhat higher densities in the bulk probably because 
the density itself does not depend on the interaction in 
a uniform Bose gas. Recent experiments have ob- 
tained stable condensates at densities corresponding to 
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Thus while the Thomas Fermi density nTF(O) at the cen- 
ter of the trap is accurate, the TF energy is a poor ap- 
proximation in this density regime. This is because the 
TF approximation underestimates the density at large r 
and the high energy particles at large r contribute signif- 
icantly to the energy. 

In Fig. 2 we see that the present Monte Carlo E/N lies 
above the Gross-Pitaevskii i;/7V, by 2.5% at TV = 20, 000. 
A difference could arise for two reasons. 

Firstly, the MC energy is the E/N of the whole Bose 
gas while the GP equation calculates the E/N of the con- 
densate only. To investigate the effects of depletion on 
E/N we note, comparing (28) and (30), that the frac- 
tional change in energy arising from depletion, SE/E, 
can be related to the fraction of atoms excited out of the 
condensate as, 

SE/E = (128/15)(naV7r)^ = il6/5)6N/N (33) 



FIG. 10: Energy difference between present VMC and GP 
results, SE/N , as a function of Na/a^o for four different values 
of a/ant — 1, 5, 10, and 20. The lines are fits of SE/N = 
m{a / aho)'^ N + 6 to the data points. 



« .85 and effects resulting from depletion are expected 
to be significant. 

For a uniform Bose gas, the MF, Bogoliubov expression 
for the energy including the leading correction arising 
from depletion is [D 



E ^ 3. 128, na^ 
- 1 + — 

iV 15 TT 



(30) 



The corresponding TF expression for bosons in a spher- 
ical trap (A = 1) is 



E 5 

— — /iTF 1 

N 7^ ^ 



77r^ nTF(0)Q^ ^ij 



(31) 



where fixF = 5?iwj^(15iVa/a/io) 5 and from ( p6| ) and ( |27| ) 

nTF{0)a^ = {15i/8TT)Ni{a/aho)^ 

The above expressions may be used to estimate the 
energy of bosons in an anisotropic trap by replacing a^o 
with the geometric mean Ug ~ (02,0^02)3 which is ag — 
aho^~^ for the trap discussed in Fig. 2. The energy of 
independent, non-interacting bosons in this anisotropic 
trap is 

(E/N) -> Eho = huJhoil + A/2) = hcjho{2AU) (32) 

{ujho = oJi_)- In Fig. 2, [E/N) has the non-interacting 
boson limit at — > 0. As N increases, {E/N) increases. 
As seen from (31), the mean field energy in the TF ap- 
proximation (E/N)tf = 5/iTi?/7 increases with N as 

{E/N)tf oc {Na/aho)i. The {E/N)mc in Fig. 2 follows 
the dependence approximately. However, direct evalua- 
tion of (31) shows that {E/N)tf underestimates the en- 
ergy significantly, by 25% at N=20,000 (Na/0,,0 = 86.6). 



in the bulk. The corresponding equation for the trap is, 
from (29) and (31), 



6E/E = {7/5)6N/N. 



(34) 



These are lowest order expressions. In Table 2 we list 
the depletion of the condensate SN/N for the bosons 
in the anisotropic trap considered in Fig. 2 for N = 
5000,10,000 and 20,000 predicted by (||). We expect 
these predictions to be accurate since nxFiO) is reliable 
and nTF{0)a^ is small. The interaction depletes the con- 
densate 0.5% at = 20, 000. From the above connection 
between SE and SN, the energy {E/N) including deple- 
tion is expected to he 0.7-1.6% above the energy of the 
condensate. On this basis, VMC and GP energies are 
consistent. 

To explain the connection between the difference in 
VMC and GP energies and depletion of the condensate 
more fully, wc have plotted S{E/N) — {E]\{c ~ Eqp)/N 
vs Na/aho in Fig. 10. From (31), the difference from the 
mean field energy arising from depletion is 



S{E/N)tf = {5TTiiTF/8)inTF{0)a^ /tt)'- 
oc N's{a/aho)'^ 



(35) 



N 


Na/aho 


nTF(Q)a^^ 


[ — )MF 


( — )MC 


5000 


21.6 


1.1 X 10"^ 


0.37% 


1.1% 


10000 


43.3 


1.5 X 10"'' 


0.43% 


1.8% 


20000 


86.6 


2.0 X 10"^ 


0.50% 


2.5% 



TABLE II: Depletion of the condensate {5N/N)mf calculated 
using Mean Field expansions ( |2£| ) for the Bose gas considered 
in Fig. 2 for N= 5000,10,000, and 20,000. The {5E/E)mc 
is the fractional difference between the present MC energy 
of the gas and the GP energy of the condensate, SE = 
Emc — Egp- This difference is consistent with the depletion 
of the condensate and mean field expressions of SE/E arising 
from depletion. 
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in the TF limit. In Fig. 10, we show the difference in 
the VMC and GP energies, S{E/N), verses the product 
Na/afio- This difference is obtained from Fig. 4 for four 
values of a/uho — 1,5,10,20. The lines in Fig. 10 are 
fits of 6{E/N) — TUN'S [a/ aho)^ + 6 to the data where 
the fitting parameters m and h where allowed to change 
for each a/aho value. The good fit of these lines shows 
that 5{E/N) reflects the dependence on N expected for 
a difference in energy arising from exciting bosons out of 
the condensate. Thus the difference in {E/N)mc from 
{E/N)gp is consistent in magnitude and dependence on 
N and a/uho with that expected for an {E/N)mc that 
includes bosons both in and above the condensate while 
{E/N)gp is the energy of bosons in the condensate only. 
Thus we believe the MC and GP energies differ chiefly 
because the MC includes "excited" particles while the 
GP energy does not. 

Secondly, the present VMC E/N is a genuine up- 
per bound for the whole gas and could lie above the 
whole gas energy. We have not tested the sensitiv- 
ity of E/N to different choices of the trial wave func- 
tion. In addition, while the present pair Jastrow func- 
tion (||) is exact in the dilute limit it does not contain 
any variational parameters and is therefor not optimized 
for trapped hard spheres at higher densities. As a re- 
sult, at least some of the difference in {E/N)mc from 
{E/N)qp could arise from the present choice of trial 
wave function. For example, Fabrocini and Polls (FP) 
have evaluated the whole E/N for bosons in a spherical 
trap — ay = az = 4.33 x 10~'^a/io using correlated ba- 
sis function and hypernetted-chain (HNC) methods [^. 
Both methods provide estimates oi E/N which lie be- 
low the present VMC E/N. At N = 10^ the density 



(n(O)a^ = 2.5 x 10 ^) in this trap is similar to that for 
N — 2 X 10* in the present elliptical trap. At this den- 
sity the HNC E/N lies 0.8% above the GP energy of the 
condensate while the present VMC E/N is 2.5% above 
GP. FP also consider a mean field model incorporating 
quantum depletion which predicts an increase in E/N of 
1.2% above the GP result at this density. Generally, the 
HNC energy in the trap lies below the energy expected 
for depletion and may be too low. For example, the HNC 
energy for a uniform gas of bosons lies above the energy 
expected for depletion. Definite resolution of these dif- 
ferences awaits a model independent evaluation oi E/N 
by Diffusion Monte Carlo methods. 

Finally, an important result of the present MC eval- 
uation is that as the boson density na^ increases, the 
condensate gradually moves from the center of the trap 
to the edges of the trap as shown in Fig. 6. At large 
na^, the condensate is at the edges of the trap. In this 
limit, the depletion of the condensate is large and the 
condensate seeks the regions of lowest total density which 
are at the surface of the trap. Both the condensed and 
uncondensed atoms must be included in the calculation 
to obtain this effect. This result is consistent with the 
calculations in liquid '^He droplets |l^ which find the 
condensate concentrated at the surface of the droplet. 
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